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A GENERAL METHOD FOR CALCULATING THREE-DIMENSIONAL COMPRESSIBLE 
LAMINAR AND TURBULENT BOUNDARY LAYERS ON ARBITRARY WINGS 


by 

Tuncer Cebeci, Kalle Kaups, and Judy A. Ramsey 
McDonnell Douglas Corporation 

SUMMARY 

A very general method for calculating three-dimensional compressible 
laminar and turbulent boundary layers on arbitrary wings is described. 

The method utilizes a nonorthogonal coordinate system for the boundary- 
layer calculations and includes a geometry program that represents the 
wing analytically, and a velocity program that computes the external velocity 
components from a given experimental pressure distribution when the external 
velocity distribution is not computed theoretically. The boundary-layer 
method is general, however, and can also be used for an external velocity 
distribution computed theoretically. 

The boundary-layer method accounts for all the geometric parameters of 
the coordinate system. The Reynolds shear-stress terms are modeled by an 
eddy-viscosity formulation developed by Cebeci. The governing equations are 
solved by a two-point finite-difference method used earlier by Keller 
and Cebeci for two-dimensional flows and later by Cebeci for three- 
dimensional flows. 

o 

Several test cases are computed by this method and the results are 
checked with other numerical calculations and with experiment when available. 
A typical computation time (CPU) on an IBM 370/165 computer for one surface 
of a wing v/hich roughly consist of 30 spanwise stations and 25 streamwise 
stations, with 30 points across the boundary layer is less than 30 seconds 
for an incompressible flow and is a little over this for a compressible 
flow. 
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SYMBOLS 


A 

b 


c 

c ^c 


Cf 


n 


C 

c p 

E 

f 

f' 

g 

g' 

G 

H 


i t j , k 


Van Driest damping length 

C(1 + e m ) or a constant for scaling the spanwise coordinate y 
P e /p or local chord length 

local skin-friction coefficient based on the Cartesian component 
of wall shear stress vector in x-coordinate direction 

local skin-friction coefficient based on the Cartesian component 
of wall shear stress vector normal to the x-y-plane. 

PU/(P e P e ) 

p 

pressure coefficient, 2(p — P^/Cp^uJ 
total enthalpy ratio, H/H e 
transformed vector potential for ip 
u/u e 

transformed vector potential for <j> 

v, / u ref 

E* 

total enthalpy 

unit vectors in x,y,z-directions of the Cartesian coordinate 
system in which the wing is defined. 


Kl=K gi , 

K 12’ K 21 

t 


M. 


n 


geodesic curvatures 
geometric parameters 

curvature vector of (3-D) coordinate line 
modified mixing-length or a scaling constant 
parameters in transformed differential equations 
free-stream Mach number 
unit outside normal vector to the surface 
static pressure 
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Prandtl number 


total pressure 

(x,y,z) position vector for point on surface 

u e s -|/v e# Reynolds number 

arc length along x-coordinate line 

time 

unit tangent vector along (3-D) coordinate line 
component of velocity vector in x-coordinate direction 
reference velocity 

resultant velocity at the edge of boundary layer 

total or resultant velocity 

derivative 3u/3x 

friction velocity 

free-stream velocity 

component of velocity vector normal to the surface 
component of velocity vector in z-coordinate direction 

3V//3Z 

5 Reynolds stresses 

independent variable in chordwise direction, x = 4> — 4> 0 
(attachment line) 

independent variable normal to the surface, is equal to the 
normal distance. 

independent variable in spanwise direction, z = y/b 

Cartesian coordinate system used for wing definition 

local geometric angle of attack of wing section chord lines 
with respect to x-axis (a(y) expresses wing twist) 

constant in outer eddy viscosity equal to 0.0168 



p 


flow deflection angle 
isentropic exponent, here y * 1.4 


s *c 


s n 


displacement thickness based on the Cartesian velocity 
components in the x-coordinate direction 

displacement thickness based on the Cartesian velocity 
components normal to the xy-pl'ane 


e (or eddy viscosity and eddy conductivity, respectively 


'H 


P' 1 9^2 »^3 
v 

% 

P 

T 

* 

Subscripts 

e 

g 

i 


airfoil ordinate in the coordinate system defined in 
figure B2* 

transformed coordinate normal to surface 

angle in tangent plane between x and z coordinate lines 

momentum thickness based on the Cartesian velocity components 
in the xrcoordinate direction 

momentum thickness based on the Cartesian velocity components 
normal to the xy-plane 

local sweep angle, measured between normal plane to freestream 
velocity vector and z-coordinate line 

dynamic viscosity 

parameters in transformed energy equation 
kinematic viscosity 

airfoil abscissa in the coordinate system defined in figure A2 

density 

shear stress 


stretching variable defined by eq.(B5) (x-coordinate in 
differential equation is x = <f> — <j> Q (attachment line)) 

two-component vector potential, eq. (33) 


outer edge, effective 
geodesic 

at chordwise input stations 
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i >o 

j 

£ 


inner and outer regions for eddy viscosity 
at spanwise input stations 
leading edge 

p chordv/ise stations at which boundary layer is computed 

r wi ng root 

s spanwise stations at which boundary layer is computed 

t wing tip or "total" 

w wal 1 

® free-stream conditions 

bars denote Cartesian coordinate system 

primes denote- differentiation with respect to n 
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INTRODUCTION 


In recent years several methods have been developed to compute three- 
dimensional laminar and turbulent boundary layers. For example, in [1], 

Hunt, Bushnell and Beckwith developed a finite-difference method for analyz- 
ing compressible turbulent boundary layers on a blunt swept slab with leading 
edge blowing. In that reference, the Reynolds shear stress terms were 
modeled by using a mixing-length expression which was based on a generaliza- 
tion of the two- dimensional model used earlier by Bushnell and Beckwith [2], 
This model has been adapted by Adams and was used to compute several semi- 
three-dimensional turbulent boundary layers (see, for example, [3]). In [4] 
Bradshaw extended his two-dimensional method to compute three-dimensional flows 
past infinite swept wings and obtained good agreement with experiment. His 
method differs from the previously mentioned ones in that his modeling of the 
Reynolds shear-stress terms does not use the mixing-length concepts; rather, 
it uses the turbulence kinetic energy equation, a model that has been used 
extensively by Nash and his associates (see, for example, [5], [6] ). 

Recently Harris and Morris [7], Kendall et al. [8], and Wortman [9]* 
developed methods to compute full three-dimensional compressible laminar 
and turbulent boundary layers. Their methods also use mixing-length, 
eddy- viscosity concepts to model the Reynolds stresses. Except for the 
methods of Kendall et al . and Wortman, none of the above-mentioned methods 
have been applied to flow problems over real aerodynamic configurations 
such as wings, and the problems of calculating three-dimensional laminar , 
and turbulent boundary layers on such bodies have never been investigated. 

In the present report, we discuss a general method applicable to 
three-dimensional compressible laminar and turbulent boundary-layer flows 
on arbitrary wings. The method uses an eddy-viscosity formulation devel- 
oped by Cebeci [10] and uses a two-point finite-difference method developed 
by Keller and Cebeci [11] • So far this method has been tested for most of 
the available experimental data on various three-dimensional flows and has 
been found to give accurate results (see, for example, [10], [12], [13]). 



The three-dimensional method presented here has a number of desirable 
features. First, it uses a very convenient coordinate system to do the 
boundary-layer calculations. This is a nonorthogonal system that eliminates 
most of the computational difficulties associated with the orthogonal systems 
used for wing problems. To illustrate this significant feature, let us con- 
sider an orthogonal system in which the orthogonal s are constructed at 
constant percent-chord stations. With this system, as shown in figure -1, 
the orthogonal s started from the wing-root congregate at the nose, leaving 
large portions near the trailing edge uncovered. This is especially true 
for a wing with a sharp trailing edge. A round trailing edge rectifies the 
situation somewhat but there is still a large area of the wing where the 
orthogonal s are sparce. 

The coordinate system in figure 1 was constructed with the polar angle 
$ - x at the root section as the other surface coordinate (see figure B2). 

As is seen from figure 1, there are computational difficulties at the trail- 
ing edge. To show this, consider fi-gure 2 in which the surface coordinate 
network x and z is obtained by extending the surface coverage with the 



Figure 1. An orthogonal system for the wing. 


2 


A" 



Figure 2. Wing in the x and z plane. 


dashed lines. Here AA" is the stagnation line, AB the root section, and 
D is a point on the trailing edge. Starting from the initial lines, the 
boundary layer can be calculated along the line BC" including the root. 
However, the point D cannot be obtained in a straightforward manner. This 
is also true for the rest of the trailing edge points D ' , D" and D WI . 

Another possible coordinate system can be obtained by representing the 
wing by one or more separate conical surfaces. Figure 3 shows such a 
representation. Here the wing panels ABDC and CDFE form two conical 
surfaces with apexes at P and Q, respectively. The shape of the panel 
ABDC and the coordinate system in the developed plane are shown in figure 4. 
The initial lines are AC and AB. AC is the stagnation line and AB is 
the wing fuselage junction. Calculations can be started at corner A. A lin 
ear coordinate transformation can be used to avoid marching into the negative 
r direction. Such a coordinate system without taking the thickness into 
account ( this amounts to representing wing sections by flat plates ) was 
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Figure 3. Representation of the wing by conical sections. 



Figure 4. Notation for conical section ABDC and the coordinate system in the 
developed plane. 
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used by Nash and Scruggs [6], The disadvantage of this coordinate system 
is the difficulty of doing calculations in the overlap region. This is a 
major difficulty and with the assumption for zero thickness, it is a poor 
and inaccurate coordinate system for boundary-layer computations. 

A second major desirable feature of our method is the geometry subprogram 
which is employed to calculate the coordinate system and its geometrical param- 
eters. By means of this program, since we are using a body-oriented coordinate 
system, we represent the wing analytically and calculate the geodesic curva- 
tures and the metric coefficients once and for all. In this way we eliminate 
the need for calculating the coordinate system, say, with each angle of attack, 
as in the case with a streamline coordinate system, 

A third desirable feature of our method is its extremely short computa- 
tion time. According to our calculations, we observe that a typical computa- 
tion time (CPU) on an IBM 370/165 computer for one surface of a wing, which 
roughly consists of 30 spanwise stations and 25 chordwise stations with 30 
grid points across the layer, is less than 30 seconds for an incompressible 
flow and slightly more for a compressible flow. There are several reasons 
for such a small computation time. First, we use transformed variables that 
allow us to take a few spanwise and chordwise stations. Second, our coordi- 
nate system is an '’optimum’' coordinate system for a wing, the flow on one 
streamwise section differs little from the other neighboring streamwise 
sections, making the numerical solutions converge very rapidly. A third, 
and maybe the most important reason, is the numerical method we use. This 
is a very efficient scheme that has been found to be very suitable for 
parabolic partial differential equations (see, for example, [14]). 

The method so far has been demonstrated for several test cases and its 
accuracy has been checked with experimental data and with other numerical cal- 
culations. The results look very encouraging. However, additional calcula- 
tions are needed to further test the method for flow conditions other than 
those studied here. 
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GOVERNING EQUATIONS 


The governing boundary- layer equations for a nonorthogonal coordinate 
system are given In references [15] and [16]. With a slight change of 
notation for compressible laminar and turbulent flows (x and z denoting 
the surface coordinates (see figure 5)), they are given by 

Continuity equation 

(puh 2 sine) + (pwh-j sine) + |y (pvh^ sine) =0 (1) 


x-Momentum equation 


p |ii + p + p V |y — p cote K-jU^ + pcsce K 2 w^ + pK-| 2 

1 2 


uw 


= _csc^ 9|£ + cotecscs |E + 3 L |£_ P irv) ( 2 ) 

TTj 3x h 2 ey \ 3y / 


z-Momentum equation 


p p pV l7- p cot9 K 2 w2 + pcsce K l“ 2 + pK 21 UW 


_ cote csce 


] 1 


<3) 


Energy equation 


u 3H, w 3H . ~rr 9H _ 3 

p Ef33r +p T^37 pV ^y"57 


U 

Fr 


ly + - Ff) ly ( t ) - pv H 


(4) 


Here pv = pv + pV and h-j and h 2 are metric coefficients. The latter 
are functions of x and z, that is. 


h ] = h 1 (x,z) , h 2 = h 2 (x,z) (5) 

Also, e represents the angle between the coordinate lines x and z. For 
an orthogonal system 0 = tt/ 2. The parameters K-| and K 2 are known as 
the geodesic curvatures of the curves z = const, and x = const., respec- 
tively. They are given by 
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K = 1 

1 h-jh 2 sine 


[h (h 2 cose > ~ Sr]’ K 2 ■ H pr -' sine " [k < h l cose > " T, r] 

( 6 ) 


ah. 


The parameters K-j 2 and K ?1 are defined by 


21 


K 


1 


12 sine 


h + b[w) + cose ( k 2 + htIf ) 


(7a) 


K 21 ' iW ~ ( K 2 + l^ff) + cose ( K 1 + Ej- If) 


(7b) 


u t represents the total velocity within the boundary layer and is given by 


u^ = (u^ + w^ + 2uw cose) 


1/2 


(8) 


At the edge of the boundary layer, (2) and (3) reduce to 


u e 9u e . w e 9u e 


p e l TTj" alT + sT ~ cot0 K l u e + csce K 2 w e + K 12 u e w e 


_ _ CSC 9 3 p cote CSC B 3 p 
3 X 32 

/u e 3W e W e 3W e 2 2 \ 

p e yFj" ax - + h^ Tz cote K 2 w e + csce K l u e + K 21 u e w e ) 

2 

_ cote csce 3p esc e ap 


'1 


3X 


h 2 3Z 


( 9 ) 


(10) 


The boundary conditions which we shall consider for (1) - (4) are: 



y * <5 u = u e (x,z), w = w e (x,z), H = H e . (lib) 

The solution of the system given by (1) - (4) subject to (11) requires 
closure assumptions for the Reynolds stresses, -pu' v 1 , -pw ( v' and -pv'H T * 
They also require initial conditions on two intersecting planes. Here we 
consider several options of starting initial conditions on these planes. When 
the solutions start at the leading edge of the wing, we use the stagnation- 
line equations. These equations are obtained in the present coordinate 
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system be examining the velocity components in the coordinate system shown 
in figure 6 where we see that 

U = u sine 
w = w + u cose 

At the stagnation line u = 0. Therefore, u = 0 and w = w. Thus, equations 
(1), (2), and (3) reduce to the stagnation-line equations if we set u = 3p/ax = 0. 
Then the x-momentum equation becomes singular on that line. However, differenti- 
ation with respect to x will yield a nonsingular equation. After performing 
the necessary differentiation for the x-momentum equation and taking advantage 
of approximate symmetry conditions (aw/ax = av/ax = a u/ax = 0) and using (9) 
and (10), we can write the governing stagnation-line equations as: 



Figure 6. Sketch showing the coordinate velocity components in orthogonal 
and nonorthogonal systems. 


Conti nui ty: 


ph 2 sine u x + (pwh-j sine) + |y (pVh-j h 2 sine) = 0 
x-Momentum: 

2 


(12) 


3U 
W X 


3U. 


u Vtt w Q au va 
xe . e xe 


TTj~ ' + p KJ’ az” + pv ay~ + pK 12 wu x = p e \~K^ + h^ “ + K 12 w e u x e/ 


9y 


au 

u ay*- - p(u T v r ), 


(13) 


I 


z-Momentum: 


w aw 


aw 


2 _ 


r w rt aw rt 
e e 


p TiTaz + pv W p cot0 ^ = p e Uraz-- cote K 2 w e 


W ( v " p " r * r ) 


Energy 

w aH . — aH _ a 
p 9z pv ay "ay 


u aH . /•. 1 \ 

Pr liy + * \ — Fr/ 


. (A' 

w\r. 


— pv'H 


( 14 ) 


(15) 


Here u„ = 3u/3x, u v „ * 3U../3X and total velocity u + = w. These equations 

X X0 0 t 

are subject to the following boundary conditions 


II 

o 

u x , V, w = 0, 

(3H/3y) w - 0 

(16a) 

y = <s 

U x - u xe (x.z). 

w = w g (x,z), H * H fi 

(16b) 


When the solutions start at the root or the tiD of the wing, excluding 
the leading edge, we use either the chordwise attachment-line equations 
(plane of symmetry) or the infinite-swept-wing equations. Note that these 
equations are just approximations to the governing equations in those two 
regions. The attachment-line equations are obtained in a manner similar 
to the stagnation-line equations. On this line w and 9p/9z are zero 
making the z-momentum equation singular. However, differentiation with 
respect to z yields a nonsingular equation. After performing the necessary 
differentiation for the z-momentum and taking advantage of the approximate 
symmetry conditions (9u/9z = 9v/9z = 9 w/gz = 0) and using (9) and (10), 
we can write the governing chordwise attachment-line equations as: 


Continuity: 


9 X 

x-Momentum 
U 9U . — rr 9U 


r\ _____ 

(puhg sine) + ph-| sine w z + ^ (pvhi hg sine) = 0 (17) 


2 . I u e 9u e 


p Ef 3x + pV 3y - p cote K 1 u = p e ytvj" 3x cote 

+ ly ( v ly - P ‘ 7V ’) ^8) 
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z-Momentum 
3W 


3W ? 3W /U 3W W > 

TTj” 3x + w z + p v 3y + p ^21 UW 2 ” p e ( Ej~ 3x + h 2 + *^21 u e w ze 


a_[ 

ay 3y 


pfw 7 ^. 


Energy 


u 3H . 3H _ 3 
p Ti^ 3X + pv ay "ay 


/ d. \ 

u 3H . /, 1 \ a f t \ -nrr 

Fr W y V 1 “T'rj ?y \2“/ _pv H 


(19) 


( 20 ) 


Here w z = 3w/az, w ze = 3w e /az and total velocity u t =.u. These equations 
are subject to the following boundary conditions: 


y = 0 

y = 6 


u = v = 0 


w z = 0 


(aH/ay) w = 0 (21a) 


u = u e (x,z) w z = w ze 


H = H. 


(21b) 


The infinite swept wing equations are obtained by neglecting the spanwise 
variation of u, v, w and H. They are given by: 

Continuity 


§x (puh 2 sine) + |- (pvB-jh 2 sine) = 0 


( 22 ) 


x-Momentum 

p F“fx + ^ fy ~ p cote K i u2 + p csce K 2 w2 + P K -j 2 uw 


gy p uui/ v i\-j u 1 yj i\2” ■ pi\-| 

(§■ W- ~ cote K l u e + csce K 2 w e + K 12Ve) + ly ( p J? ~ p “’* r ) 


= P, 


z-Momentum 

p F _ lx' +p ^ly _p cote K 2 w2 + p csce K l u2 + pK 21 uw 

/ u e 8w e 2 2 \ 

= p e( H7 air - cote K 2 w e + csc9 K l u e + K 21 u e w e) + 


(23) 


— ( 
ay V 


aw 


w 3y- PW'V 1 


(24) 
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Energy 


u 9H 
p TTj~ ax 


+ 



3 _ 

3y 


9H . (i 




( 25 ) 


The boundary conditions are the same as those given by (11) except that u e 
and w e are independent of z. 
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CLOSURE ASSUMPTIONS FOR THE REYNOLDS STRESSES 

For turbulent flows, it is necessary to make closure assumptions for 
the Reynolds stresses, -pu T v r , -pw'v* and -pv'H' . In the study reported 
here we satisfy the closure conditions by using the eddy-diffusi vity concept 
and relate the Reynolds stresses to the mean velocity and total enthalpy 
profiles by 


. 3U 

" pu , v = pe m ay 

9 

,■ i 3W 

- pw v = P£ m sy 

(26a) 

— nrr — a H 

-pv H - PS H jy 



(26b) 


By means of "turbulent" Prandtl number, Pr t (= e m /e H )» we can also write 
(26b) as 


i ,, i m 

-pv H - p p^r~ 


3H 

ay 


(26c) 


We use the eddy-viscosity formulation of reference [10], and define e m by 

two separate formulas. In the so-called inner region of the boundary layer, 

(e ) is defined by the following formula 
v nr 



(27) 


where 


L = 0.4y[l - exp(-y/A)] 



In the outer region e m is defined by the following formula 


(28a) 

(28b) 

(28c) 


(fij 

m . 


0.0168 



- u t )dy 


(29) 
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where 


= (Ug + Wg + 2UgWg coso) 1/2 

(30a) 

= (u 2 + w 2 + 2uw cose) 1/2 

(30b) 


The inner and outer regions are established by the continuity of the eddy- 
viscosity formula. 


According to a number of studies, the turbulent Prandtl number is 
relatively constant across the boundary layer (see reference [14], for example). 
In our study, we shall take it equal to 0.90. 
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TRANSFORMATION OF THE GOVERNING EQUATIONS 


The boundary-layer equations can be solved when they are expressed either 
In physical coordinates or In transformed coordinates. Each coordinate has its 
own advantages. In three-dimensional flows, where the computer storage and 
time becomes quite important, the choice of using transformed coordinates 
becomes necessary as well as convenient because the transformed coordinates 
allow large steps to be taken in the streamwise and spanwise directions. In 
addition, they remove the singularity the equations have in physical coordinates 
at x s 0 and z = 0. 

(4). 
(31) 


(32) 


(33) 

Using these transformations and the relations given by (9), (10) and (26a, c), 
after some rearranging, we get 

x-Momentum 

(bf” ) + m-|ff" — m 2 (f ') 2 — mgf'g' + mgf'g - m 8 (g ‘ ) 2 + m-| -jC 

= *10 f If - f " w) + m 7 (s’ If - f ” H)< 34 > 


Let us consider the transformation of the equations given by (1) 
In this case we first define the transformed coordinates by 


x = x, z = z, 


U e \1/2 

tin ■ --T- 1 pdy. 


v p e u e s l 

and introduce a two-component vector potential such that 
puh 2 sine - 1^- , pwh-| sine = |y 

^h 2 sin e -(fftff) 

In addition, we define dimensionless ip and cj> by 

1 /2 

ip = (p e y e u e s-|) h 2 sine f(x,z,n) 


i = / h i dx 

0 


V2 


= (p e v» e u e s 1 ) u ref /u e h-j Sineg(x,z,n) 
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z-Momentum 

(bg") + m^fg" -n^f'g' -m 3 (g') 2 + m 6 gg" -m 9 (f') 2 + m 12 c 


= mn 


( f - 


10 I' 3X 


n» 'if 


ax 


) + m 7 ( 


g' <L9l_gl' M 

M , 3Z .5* DZ 


■K 


Energy 


{ Wl E')’+ p z E*. (f H-E' f |)+* 7 (g- If- E' ff) 

Here primes denote differentiation with respect to n and 


r ■ u/« e . 


9' - w/u ref . 


b = C(1 + e*}, C = , C = — - 

v m d_u_ p 


fVe 


E = H/H 


e m - e m /v 
m m 


The coefficients m-| to m-j 2 , p-j to u 3 are given by 

i—l/2 


1 / s i ^ u e\ s l( p p p p) % 

= 7y + Jy^Sr) + sine ~ ax" (h 2 sin9 


ITU = 


'1 


*■ 1 e 

] 3 = ” S 1 


rn. = S-, 1 C 

4 I £ 


m 10 = TTj” 


m ll. S 1 



u e 3X } 

3U e 

“ S 1 K 1 ' 

ax 

cote 

u ref 

K 2"V 


S' 

1 » 

m 5 = TT 

i 


u e s l h 2 

ref 

* m 8 : 

u e 

9 


1 

9U 

. — _ + 


2 u‘ 


3u e + 

„ u ref 

al“ + 

K 12 s 1 u 

e 

-1 3 

/ 

1 az 

fWi 


/u ref\ 2 . 


(35) 

(36) 

(37a) 

(37b) 


(38) 


sine 


'ref 


S 1 u ref / u ref\ 2 u e 

m 7 s hJU ~ 9 m 8 ” s 1 K 2 csc0 \ %~) 9 m 9 = S 1 K 1 CSCe u pef 


e n l 


w au /w \ W P 

-T- sz^-cote R, + csce Uf) + K ^ 
Ugh 2 \ e/ ' e 
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m 


S 1 / u p aw e w e 9w e 2 2 , , t 

12 - I nr— ( h7 ST + TCT zT - cote K 2 w e + csc8 K 1 u e + K 21 w e u e ) < 38 > 


e ref V'l 


C /t. + Pr \ . , 

y l Pr V e m Pr t /’ y 2 m l f m 6 9, 

2 . 

"3 = C IT 0 - Pr) f,f " + g'9" (-U^f + “u^ cose ' 9" + g'f") 


(39) 


To transform the governing stagnation-line equations, we define the 


transformed coordinates by 


x = x. 


2 = z. 


dn = 


u xe \’/2 


Pe^e 11 ! 


pdy 


(40) 


the two-component vector potential by 

pu x^2 S1n0 = fy » pwh-j sine = fy 


pvh-jh 2 sine = + |jj 


and the dimensionless ip and <j> by 


(41) 


(42) 


1 / 2 

^ = ^ p e y e u xe h l^ h 2 sine f ( z » p ) 

4> = ( p e y e u xe h l^ /2h l sin0 u^~ 

With these new variables, for laminar flows (b = C), the stagnation-line 
equations can be written as 

x-Momentum 

(bf 11 ) + ff” - (f 1 ) 2 -mgf'g' + m 6 f"g + m^c = m y (g* || f" ||-j‘ (43) 


z -Momentum 


(bg n ) + fg" ~m 3 (g‘) 2 + m 6 gg" + m 12 c = m y (g‘ — g" ff-) 


(44) 
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Energy 


(uiE-)'+ » 2 E* + p| - m 7 (g* |f ~ E* ff) 


Here mg* mg, mg, m-j-j and are * 


m 3 = “h-jKg cote u 


ref 


xe 

u u ref /„ . ^2 3u xe 

m5 1 vT V 12 ^ 3Z 


m, 


! 6 ' [<Pe l 'e u xe h l) 1/2 h 2 sin0 ^' fz 


,-1 3 


VPe^e^e^ h l s1ne u 


ref 


h l u ref 


w_ 


1 . h 1 w e ^xe , . "e 

mi1 1 V 12 

xe 


h l w e 
m 12 = ‘Hr Ti~u 


ew g h-j cote I^Wg 


2 u xe u ref 92 u xe u ref 


xe J 


(45) 


(46) 


The coefficients y-j » P2» Ug are: 


y l ‘ Pr 


y 2 


= f + 


ni 6 g 


_ c 

y 3 H 


(l -^)u 2 ef g'g"(47) 


To transform the governing infinite-swept-wing equations, we use the 
transformation given by (31) to (33) and write (23) to (25) as: 


x-Momentum 


(bf")' + m, ff" -m 2 (f') 2 - m 5 f 'g' — m 8 (9 ' ) 2 + m 7) c 

. S 1 l f , af ' f „ H \ 

hj \ 3 X 3 X / 


(48) 
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z -Momentum 


(bg") + m-jfg" -n^f'g' -m 3 (g') 2 -m g (f') 2 + m 12 c 




'laL — *1 

3X a 3X 


( 49 ) 


The energy equation is the same as the energy equation for the general case, 
(36), except that the right-side of (36) is: 


- S 1 lf> 21 r. 3f\ 

■ y E lx) 


(50) 


The coefficients and y 3 remain unchanged. However, u 2 now becomes 

equal to m-| f . 


The definitions of the coefficients m-j to m-| 2 in (48) and (49) are 
the same as those in (38) except now 

_ i/ u ref 
m 5 " K 1 2 S 1 1T~ 


m 6 = 0 


m ll S 1 


1 3u e 

u^fTf 3X- - cot9 K 1 + csce K 2 


m 2 + k ^ 

u e/ 12 u e 


(51) 


12 u e u ref V'l 


u e 3w e 2 2 

H7 3x cote K 2 w e + csce K l u e + K 21 w e u e 


To transform the governing chordwise attachment-line equations we use 
the transformed coordinates given by (31) and define the two-component vector 
potential by 


p uh 2 sine 


3y 


w z ph-j sine 


B4> 

3y 


pvh-|h 2 sine 



+ <f> 


) 


(52) 


with ip and <f> still given by (33). With these variables, and with the 
relations given by (26a,c), the chordwise att ac hment-line equations can be 
written as 
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x- Momentum 


(br)' + m-j f f " — m 2 (f 1 ) 2 + m 6 f"g + m^c 

_ /ri 3f 1 _ .pn 3f \ 

^ 1 3X 3 X/ 

z-Momentum 

(bg") + m-jfg" -m 4 f , g' — m 3 (g ' ) 2 + m 6 gg" -m 9 (f‘) 2 + m 12 c 

■ S ( f ' If - - 9" W) 


(53) 


(54) 


The definitions of the coefficients m-i to m-.- in (53) and (54) are the 
same as those in (38) except now 


(Tin = 


S 1 u ref 


3 " ^ u e 


m 6 = m 3 . 


m 9 = 0 

S 1 1 3u e 

m ll = PTJ' W~~ s l cote K 1 


(55a) 


m 


- i2_LL!!ze + Wze + K w 
12 u— jrlhi 3x u„ho ^21 w ze 


ref V’l 


e 2 


The energy equation is the same as the energy equation for the general case 
(36), except that the right side of (36) is the same as (50); the coefficients 
vi j and U 2 remain the same but now 

.2 

-jLjf'f" (55b) 


Cu * 


The general boundary conditions for the governing three-dimensional 
boundary-layer equations are: 


T) = 0 

II 

*- 3 

It 

g w = ^ = o, 

E 1 = 0 
w 

(56a) 


f' = 1 , 

9' = w e /u ref 

E = 1 

(56b) 
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Those equations (56) apply to the general case, to the stagnation-line equa- 
tions, to the infinite-sweptrwing equations as well as to the chordwise 
attachment-line equations, except,, however,, that the edge boundary condition 
(g * a t n = Tij is different for the chordwise-attachment-1 ine equations; 

it is given by g‘ = w ze / u ref* 

■/ 

Eddy Viscosity Equations 

If we apply the transformation given by (31) to (33) to the eddy-viscosity 
formulas given by (27) to (30), we get 


. 2 

( 4 > . - (r ) Rel/2 f (°- 4 / cdn ) [ ~ ex P (- - i ) 

1 \ n A L J 


(r) 2 + (w'\ (g „)2 + 2cos6 !ref fV 


1/2 


(«£) ■ Gr) Re 1/2 (0.0168) / C 


1 + 


w\2 


/ W e\ 

+ 21—] cose 
\ u e/ 


1/2 


(57) 


( f . ) 2 + /isty (g . ) 2 + 2 !^tc 0 , 8 fV 


Here Re = n^s-j/vg and 


dr, 

(58) 


*- 


. Vl_ 0 J/4 


2F Re , 


(/ 


cdn C y 2 c 1/2 


(f") 2 + 


+ 2 cose TT f W 


1/4 


(59) 
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NUMERICAL METHOD 


We use the Box method to solve the governing equations. This is a two- 
point finite-difference method that has successfully been applied to two- 
dimensional flows by Keller and Cebeci and to three-dimensional flows by 
Cebeci. A detailed description of the method is presented in References Ql] 
and [14], For this reason only a brief description of it will be presented here. 

One of the basic ideas of the Box method is to write the governing 
system of equations in the form of a first-order system. Thus, in our case, 
the first derivatives of f, g and E with respect to n are introduced 
as new unknown functions. With the resulting first-order system and an 
arbitrary rectangular net, we use simple centered difference quotients and 
averages at the midpoints of net rectangles and net segments to get second- 
order accurate finite-difference equations. Then nonlinear difference equa- 
tions are linearized by using Newton's method and the resulting linear system 
is solved by the block-elimination method discussed by Isaacson and Keller [17], 

Numeri cal Formulation of the Momentum Equations 

In our present method we solve the two momentum equations simul taneously. 
Essentially the stagnation-line equations, infinite-swept-wing equations and 
the chordwise attachment-line equations are two-dimensional flow equations in 
the sense that these equations have two independent variables, (x, n) or 
(z,n). On the other hand, the two momentum equations (34) and (35) are three- 
dimensional flow equations for obvious reasons. The solution of the two- 
dimensional flow equations is discussed in considerable length in references 
[11], [14], [18]; for this reason we shall only discuss the solution of three- 
dimensional-flow equations, namely, (34), (35) and their boundary conditions, 
(56). 


With the introduction of new independent variables u(x,z,n), v(x,z,n), 
w(x,z,n) and t(x,z,n), the equations given by (34) and (35) can be written 
as 

f 1 = u (60a) 

u 1 = v (60b) 
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g' = w 
w' = t 


(bv) ' + m-jfv - m 2 u 2 - m g uw + m g vg - ny/ 2 + m^c = m 10 |^- - v 

*■7 (“§?->«) 

2 2 / 3 W 3 f ■ 

(bt)' + m-|ft - m^uw - m 3 w + m 6 gt - m g u + m^c = m 10 I u — ~ t 1 

*■? (»!?—* if) 


(60c) 

(60d) 


(60e) 


(60f) 


We next consider the net cube shown in Figure 7 and introduce the net 
points by 


= 0 

x n = 

Vi + 

k n 

n = 1,2, 

...» N 


= 0 

z i = 

z i-l + 

r i 

i = 1,2, 

t • • y I 

(61) 

= 0 

"J * 

j-1 + 

h j 

j - 1,2, 

• • • 9 J 




Figure 7. Net cube for the difference equations for three-dimensional flows. 
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The difference equations which are to approximate (60a) to (60d) are 
obtained by averaging about the midpoint (x n , z . , 


where, for example. 


f n , i f n , i 

_ 

FH 


“j* 1 

H j 


h j 

„n,i 

n,i 

9 j 

- 9 j-i 


h j 

-*j-5 

h j 


u 


n,i 

j-1/2 



w 


n,i 

j-1/2 


t n,i 
r j-1/2 


u 


n , i 1 

j-1/2 7 



111 '’) 

J-l 


(62 a) 


(62b) 


(62c) 


(62d) 


The difference equations which are to approximate (60e,f) are rather 
lengthy. To illustrate the difference equations to equations similar to 
(60e,f), we consider the following model equation: 


V + m i fv * m io u fj + V fr 

The difference equations for this equation are 


Vj - V.. 


u„ — u 


j-l , /_ \ n-1 /2 rzr\ _ ;n-l/2 - " / u n' u n-l 

TT + ("l^ 1-1/2 j-1/2 ^ m 10^i-l/2 u j-‘ 


+ ^ m 7^i-l/2 w j-l/2 


■1/2 v " , 10 / i-l/2 j-1/2 \ k p 

i n "l/2- ^ U i ~ U i-t 

where, for example, 

v, -i (v "- 1 + v "* 1 - 1 + vr 1 * 1-1 + vr 1,1 ) 


'j • T ' v j 

5 n " T (“j*' + -S’ 1 ' 1 + + ^l* 1 ) 


(63) 


(64) 
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«,-i 

The boundary conditions for. the system of equations given by (56), ^t x = x 

n,i 


and at z = are: 


rn,i ^ 


- 0, gg’ 1 = 0, u^“ = 0, wH*' = 0, u^ 1 ’ 1 =1, w; 


,n»i = 


n,i _ 


,n,i . 



(65) 


If we assume u^ 1 ’ 1 ' 1 , v"' 1 - 1 ' 1 , g^ 1 - 1 " 1 . wj- 1 * 1 ' 1 . tj- 1 * 1 ' 1 ). 

CfJ* 1 - 1 . ii!}* 1 ' 1 . gj’ 1 " 1 , w"* 1 - 1 , t"* 1 '- 1 ), and (fj' 1 * 1 , u^.\ 

v j ' ,1 » g”" 1 ,n , to be known for 0 <_ j then the differ- 

ence equations (62) and the difference equations for (60e,f) along with (65) 
yield an implicit nonlinear algebraic system of 6J + 6 equations in as many 
unknowns (f 1 ?, u!J, v 1 ?, g 1 ?, w 1 ?, t^). We solve this nonlinear system by means 
of Newton’s method. The resulting linearized system is then solved by 
using the block elimination method discussed by Isaacson and Keller [17]. 


Numerical Formulation of the Energy Equation 

The numerical formulation of the energy equation Is very similar to the 
formulation described for the momentum equations. To reduce (36) to a first- 
order system, we introduce a new independent variable G(x,z,n) and write 
(36) as 

E ' - G (66a) 

( Wl G) ' + y 2 G + = m 1Q (u || - G + m 7 (w §§ - G §§ ) (66b) 

The difference equation -for (66) is written again by averaging about the 
midpoint (x n , z^ , similar to those given by (62). The difference 

equation for (66b) is written similar to (64). The boundary conditions in 
(56) for the energy equation are: 


?5 



(67) 



p n»i 


= 1 


The resulting algebraic system of 2J = 2 equations in as many unknowns 

(E*?, g"), which is linear, is directly solved by the block elimination 
J 3 
method. 
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RESULTS AND DISCUSSION 


General Discussion of the Method and the Solution Procedure 

The present method is a very general boundary-layer method for three- 
dimensional flows. For a given pressure distribution, boundary-! ayer calcu- 
lations can be done for infinite-swept wings or finite wings. The pressure 
distribution can be either experimental or theoretical. When it is experi- 
mental, we use the approximate procedure described in Appendix A to get the 
velocity components from the experimental pressure distribution. 

The boundary-layer calculations can be started either at the stagnation 
line or at some specified percent chord away from the stagnation line. The 
initial conditions on the root section can be started by using either the 
infinite swept wing assumptions or the chordwise attachment-line equations. 
Sometimes, as shall be described later, it is necessary to prescribe the 
initial conditions at the tip section. In such a case, we use the infinite 
swept wing equations, which, as previously discussed, are an approximation 
to the flow in this region. 

The present boundary-layer method uses a body-oriented nonorthogonal 
coordinate system. A separate program is used to calculate the coordinate 
system and its geometric parameters such as the geodesic-curvature- 
parameters K-j , K^, 1C|2> l<2i and the metric coefficients h-j , hg. This 
is discussed in detail in Appendix B. 

In the present program for a given spanwise station we march in the 
streamwise direction. Since the linearized form of the equations are being 
solved, we iterate at each x-station until some convergence criterion is 
satisfied. For both laminar and turbulent flows we use the wall shear 
parameter fJJ as the convergence criterion. Laminar flow calculations are 
stopped when 
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where the value of 6-j is prescribed. A typical value of 6-j is 1CT 4 . 
For turbulent flows, the convergence criterion is slightly different; It 
is given by 


_? 

where a typical value of 6 2 is 10 . 

Whether the initial calculations start at the root section or at the tip 
section is determined by the sign of the external spanwise velocity component 
w 0 . To illustrate our marching procedure, let us assume that the wing has 
three regions defined by the sign of w g . In region 1, w 0 is positive; in 
region 2, w £ is negative; and in region 3, w g is positive (see figure 8). 
Let us also assume that region 1 starts from the leading edge. In this case, 
our calculations start at the root section. At the leading edge we use the 
stagnation-line equations for one x-station and switch to either chordwise 
attachment-line equations (if w = 0) or to the infinite swept-wing equa- 
tions. With either one of these equations we march in the streamwise 



Figure 8. Definitions of various regions on the wing for marching procedure. 
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J) 

T) 
V 'j 
> 
Of. 


direction to the net point before w e changes sign. Then we go back to the 
stagnation line and after solving the stagnation-] ine equations, we start solv- 
ing the general three-dimensional flow equations at the second x-station. This 
procedure is repeated as before for all stations up to and including the tip 
section in region 1. The calculations are continued into region 2 starting 
at the tip by solving the infinite swept-wing equations along the tip up to 
the beginning of region 3. . At that time we go back to the next spanwise 
station of region 2 and start solving the general three-dimensional flow equa- 
tions for region 2. This procedure is again repeated for region 2 up to and 
Including the root section in region 2. The calculations at the root section 
are extended to region 3 by solving either the infinite-swept-wlng equations 
or the chordwise attachment-line equations, and the procedure used for region 1 
is repeated for this region. 

It should be pointed out that our marching procedure utilizes the sign 
of w 0 . It may be argued that the marching procedure should be based on 
the local value of w. In the calculations presented here, no difficulties 
were encountered when we used the above procedure including the cases where 
local cross-flow velocity changed sign within the boundary layer. However, 
this point needs further exploration. 

The present method is developed in such a way that one can use non- 

uniform net spacings in the streanwise and in the spanwise directions. 

Across the layer one can either use a uniform grid or a variable grid 

discussed in reference [14], According to this grid, the net in the 

n-direction is a geometric progression having the property that the 

ratio of lengths of any two adjacent intervals is a constant; that is, 

h- = Kh- The distance to the j-th line is given by the following 
J i 

formula: 

nj = ( K J - 1 )/(K - 1) K > 1 (68) 

There are two parameters: h-| , the length of. the first An step, and 
K, the ratio of two successive steps. The total number of points J 
can be calculated by the following formula: 
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1n[l + (K — DCn^/h,)] 

J - TFTk — < 69 > 


For further details, see reference [14], 


In the present method the pressure-gradient parameters m-j to m-| 2 
are determined numerically from the given external velocity distributions 
u e and w e . For example, the derivative of du e /dx is obtained by using 
three-point Lagrange interpolation formulas given by (n < N) 


du a u M u 

e = fx - x ) + -S. (x 

dx n A 1 ^ x n+l x n' A 2 ^ x 


Here N refers to the last x n station and 

A 1 " K " x rt-l^ x n+l “ x n-l> 

A 2 = ^ x n “ x n-l^ x n+l “ x n^ 

A 3 = (x n + l - x n>( x n + l - x n-l } 

The derivative of du 0 /dx at the end point 


n+l _ 2x n + x n-l ^ 

< x n-Vl> 


n = N is given by 


(70) 


(71) 


du e 

where now 


TT 


N-l 


(x f 


X N-1 ^ 


( X N x N-2^ 


*7 


(2x w ~ 


x N-2 X N-1 * 
(72) 


A 1 = ^ X N-1 “ X N-2^ X N “ X N-2^ 
A 2 = ( X N-1 “ X N-2^ X N “ X N-1 ^ 
A 3 = ^ X N “ X N-1^ X N “ x N-2^ 


(73) 


Similar derivatives for du Q /dz, dw /dz, dw /dx can be written by slightly 

C “ C 

modifying the above equations. 
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Boundary-Layer Parameters 

Once the solutions are obtained, we calculate the usual boundary-layer 
parameters, which in terms of physical variables and transformed variables 
for the general case are: 

Chordwise local -skin-friction coefficient 




2x c 2 (t x + t 2 cose) 

? 2 

P co u oo Poo U co 



( 74 ) 


Spanwise local-skin-friction coefficient 


c fn 



sine 


u 


oo 



sine g" 


( 75 ) 


Cross-flow angle 



Chordwise displacement thickness 


5 


* = 
C 





I 

P 


(f » + u ref /u e cos0 9a, 5 
(1 + w e /u e COS0 J 


Spanwise displacement thickness 



( 76 ) 


( 77 ) 


( 78 ) 
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Chordwise momentum thickness e. 


oo 

-s 


pu / ! _ u_) dy = 

0 p e u e v 5 -/‘ 


(f„ + u ref /u e cose gj 
(1 + w e /u e cose) ~ 


-/ 


n2 


( fl + u re f/ u e cose 9* ) 


0 0 + w e /u e cose)' 

Spanwise momentum thickness e n 


dn 


= f / i _ 

n J . \ ,7, , 


0 p e w e 



dy = 


S 1 / u ref 
R s \ W e / 


9 »--^ / (9') 2dr > 


( 79 ) 


( 80 ) 


Here 


u^s 


R. = 


el 


S v. 


■a* 


( 81 ) 


and u/u e and w/w g are the chordwise and spanwise velocity components given 
by 


u u + w cose 


f ' + ( u re f/ u e )9' cos e 


u u e + w e Cose 1 + w e /u e COS0 


( 82 ) 


w _ w sine _ u ref w _ u ref , ,g 3 \ 

- w Q sine " w Q u^* 9 ^ ' 

w„ e . e ref e 

e 

On the plane of symmetry these boundary-layer parameters are defined by: 




Cr = <S* = * 0 

f n n n 


( 85 ) 


OO 

6* = /* /l — EH. 

C jf V p e u < 


dy --=(i„- f ) 


(86) 
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On the stagnation.! ine they are: 


Cfc = <5* = 0 C = 0 (88) 



Results for Laminar Flows 

To test our code for laminar incompressible flows, we have considered a 
test case computed earlier by Cebeci[12], by using another computer program. 
This test case consists of a three-dimensional laminar flow past a flat plate 
with attached cylinder (see figure 9). For this flow the inviscid velocity 
distribution is given by 

u e =u 4 1 + a2 ^)’ w e = - 2u =/^ W 

\ A ] / ^ 

where 

A 1 = (x - x Q ) 2 + z 2 , a 2 = -(x - x Q ) 2 + z 2 , a 3 = (x - x Q )z (93) 

Here u ro is a reference velocity, a is the cylinder radius, and x Q 
denotes the distance of the cylinder axis from the leading edge x = 0. 
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Figure 9. Flow past a flat plate with attached cylinder. 


The calculations were made for a Cartesian orthogonal system with 
u co = 3050 cm/sec, a = 6.1 cm, x Q = 45.7 cm for several x and z-stations. 
The present results with those of reference [12] are shown in Table 1. Note 
that the pressure-gradient parameters in the present method are obtained 
numerically, whereas in the method used in reference [12], they were obtained 
analytically. As can be seen from the results, there is essentially no 
difference between the earlier and the present results. All calculations 
were made with 11 points across the boundary layer and with Ax = 1.220 cm 
and az = 0.610. cm. 

In Table 1 the results under the present method represent the solutions 
obtained by the present computer program in which the pressure-gradient param- 
eters m-j to m^ are obtained numerically from the given external velocity 
distribution. The results under reference [12] represent solutions obtained 
by another computer program to which m-j to m 12 are calculated analytically. 
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Table 1. Calculated Wall Shear Parameters for the Laminar 
Post Problem. 


X 

Com) 

z 

= 0 

z = 0.610 cm 

Present 

Method 

f" 

w 

Ref. [12] 

f" 

w 

Present 

Method 

f" 

w 

Ref. [12] 

f" 

w 

0 

0.330854 

0.330854 

0.330854 

0.330854 

1.220 

0.329498 

0.329498 

0.329461 

0.329461 

2.440 

0.327973 

0.327973 

0.327711 

0.327712 

3.660 

0.326233 

0.326233 

0.325852 

0.325853 

4.880 

0.324252 

0.324251 

0.323584 

0.323588 

6.100 

0.321987 

0.321985 

0.321103 

0.321109 

7.320 

0.319416 

0.319385 

0.318078 

0.318112 


Results for Turbulent Flows 

For turbulent flows, the accuracy of the method depends on the turbulence 
model as well as on the accuracy of the numerical method. According to the 
studies reported in references [101 [12] and [13], the present eddy-viscosity 
formulation gives quite satisfactory predictions for three-dimensional flows. 

In the studies reported here, we have considered three test cases to further 
test the computer program. 

Test Case 1 

The first test case deals with the comparison of calculated and experi- 
mental results on a swept wing for an incompressible flow. The experimental 
data is due to Brebner and Wyatt [19]. However, most of this data is very 
difficult to use for comparison purposes since the location of transition is 
given at one spanwise station only. Furthermore, the relatively large probe 
used in their experiments seems to interfere with measurements at higher angles 
of attack. The "good" data is near mid-semispan of the 45-degree swept wing 
which has a 50.8 cm constant chord and an RAE 101 section with a thickness to 
chord ratio of 0.12. It has no twist and has an aspect ratio of 5. The tests 
were conducted in RAE Bedford Low Speed Wind Tunnel at a nominal free-stream 
velocity of 61 meters per second, corresponding to a streamwise chord Reynolds 
number of 2.1 x 10^. 
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For the data of Brebner and Wyatt we have made two different sets of 
calculations by using the infinite swept wing approximations for the external 
flow. With this assumption the external velocity components can be calculated 
exactly from the experimental pressure distribution. The first set of calcula- 
tions are for an orthogonal system. For the external u Q (x) distribution shown 
in figure 10 with w Q = 43 m/sec, the solutions were started as laminar at 
the stagnation point with the stagnation-line equations and were continued with 
the infinite swept-wing equations up to the transition location, which according 
to the experimental data, was specified at (x/c) = 0.35. At that location, the 
turbulent flow calculations were started and were continued up to the trailing 
edge. At the next spanwise station, the calculations were done for full 
three-dimensional flows. Overall, 30 x-stations, 15 for laminar flow and 15 for 
turbulent flows and 3 z-stations were used. The variable grid parameters were 
h-| = 0.01, K = 1.26, = 31 yielding 30 points across the boundary layer at 

the trailing edge. 

Figures 11 and 12 show the results. Figure 11 compares the experimental 
and calculated total velocity profiles at two measured x-stations, namely, 
at (x/c) = 0.80 and 0.98. Figure 12 shows a comparison of calculated and 
experimental cross-flow angle 3 , calculated from 

6 = ta,,- 1 (E)_( f _ e ) (94a) 

At the wall 6 is calculated from 

» -**"' 1 ($Mt -(*-•) i,n> 

The second set of calculations were made for a nonorthogonal system with 
hj = 0.05, K = 1.3, = 69 yielding 24 points across the boundary layer at 

the trailing edge. For the specified airfoil section, the geometric parameters 
were determined by using the procedure described in Appendix B and the velocity 
components u. and w_ were determined (see figure 13) from the experimental 
pressure distribution by using the procedure described in Appendix A. 

We note that in this case we have two regions of positive w 0 separated 
by one region of negative w_ as shown schematically in figure 8. For this 

c 
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Figure 14. Computed cross-flow profiles at different chordwise stations for the data of Brebner and 
Wyatt, nonorthogonal system, u^ = w g at x = 0. 





reason, the calculations were started at the stagnation line of the root 
section and the marching procedure described earlier was used to compute the 
flow field for three spanwise stations consisting of thirty streamwise stations. 
Region 1 is a very small region, whereas region 2 is a very large region. 

Figure 14 shows the variation of the cross-flow velocity across the boundary 
layer for various x-stations. It is interesting to note that during the march- 
ing procedure the solutions converge very fast, show no numerical problems 
and all the computed velocity profiles, when rotated to the orthogonal coord- 
inate system, agree extremely well with those results obtained for an ortho- 
gonal system. If we use bars to denote the velocity profiles in the orthogonal 
coordinate system, then the velocity profiles in the orthogonal system can be 
obtained from those in the nonorthogonal system by the following formulas 
given previously (see figure 6), that is. 


u = u sine (95a) 

w = w + u cose (95b) 


If we denote the dimensionless stream function in the nonorthogonal system by 
f and its derivatives with respect to n» then the dimensionless stream- 
function in the orthogonal system F and its derivatives are related to f 
and its derivatives by 


F' = f' 


P 1 * f 


^sin 


e s 


Similarly, 



(96) 


(97) 


Test Case 2 

The second test case is an incompressible flow past an untapered, 
untwisted, 30 degrees-swept-wing having NACA 0012 streamwise airfoil sections. 
Its semispan is 2.78 meters and its chord is 1.15meters. The inviscid velocity 
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distribution was obtained by using Hess's potential flow method [20] for an 
angle of attack of 8 degrees, giving a lift coefficient of 0.51. A total of 
16 chordwise strips on the semi span were taken and each strip was subdivided 
into 50 elements. Because the calculated inviscid surface velocities are 
given in terms of their Cartesian components, additional calculations were 
made to put these velocities into the nonorthogonal coordinate system. Calcu- 
lations were made only for the lower surface of the wing for a chord Reynolds 
number of R c = 5.88 x 10^. 

The external velocity distribution for this wing is quite interesting 
in that, except for a very small region, all the external spanwise velocity 
components w 0 are positive (see figure 15). Figure 16 shows the computed 
chordwise local skin-friction coefficients and displacement thicknesses at 
various spanwise stations. According to our marching procedure, the solu- 
tions start at the root section by using the stagnation-line equations at 
x/c = 0 and by using the infinite swept wing equations at x/c > 0. With 
transition specified at (x/c) = 0.02, the calculations continue up to the 
trailing edge since w e is positive for all chordwise stations. The pro- 
cedure is repeated for the next spanwise station, NZ = 2. The procedure 
at NZ = 3, 4, 5 and 6 are different. Since, at some chordwise stations, 
w fi becomes negative, the solutions that originate from the stagnation line 
at a given spanwise station 3 £ NZ £ 6, continue up to the last chordwise 
station where w e is positive, say (x/c) Q , and then go back to the stag- 
nation line. At NZ = 7, the procedure used for NZ = 3, 4, 5, 6 is repeated 
for the same chordwise stations. At (x/c) > (x/c) Q , however, they continue 
up to the trailing edge by using the infinite swept wing equations. At 
8 £ NZ £ 22, the procedure used at NZ = 2 is repeated. In this way all 
the chordwise stations for the spanwise stations NZ = 1 , 2, 7 to 22 are com- 
pleted. To complete the rest of the chordwise stations for NZ = 3, 4, 5 
and 6 we go back to NZ = 6 and by using the infinite swept wing equations, 
we start marching at (x/c) > (x/c) Q up to the last (x/c) station, say (x/c)-j, 
where w £ is negative. Then we go back to NZ = 5, 4, 3 and repeat the 
same procedure, but this time using the general three-dimensional flow equa- 
tions. At NZ = 3, we turn back. By using the infinite swept wing equations, 
we start marching at (x/c) > (x/c)-j. Then we go to the next spanwise stations 
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Figure 15. The wing for test case 2. The symbols denote the stations where the bound ary- layer calcula- 
tions are made. Dots correspond to stations were w fi is positive and x's correspond to 
negative w . 





Figure T6. Computed lower surface chordwise local skin-friction coefficients and displacement thickness 
at various spanwise stations for test case 2. Note that the solutions start at the stagna- 
tion line which is at x/c = 0.01, a = 1.0 (ft), a = 0.305 (meters) 
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Figure 16. 
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Concluded. 



N2 = 3, 4, and 5, arid repeat this procedure with the general three-dimensional 
flow equations, thus completing the calculation over the entire wing surface. 

j / ! 

Test Case 3 1 j 

The first two test cases considered for turbulent flows were for incompres- 
sible flows. 1 The third test case is for compressible flows. The wing selected 
for this case has a planform and dimensions similar to that given in figure 4 
of [6], except that its trailing edge is straight. It is given below as figure 
17. Also, in our system the wing root is at y = 0.76 meters corresponding to 
the z = 0 degree plane in [6]. The wing has a 5.6-meter root-chord, 1.02 meter 
tip chord; the tip is.. at y = 7.1 meters. Thickness is added by supercritical 
airfoils of 12.38-percent thickness at the root, and 9.24-percent thickness at 
the tip with parabolic variation inbetween. Twist was generated by a linear 



Figure 17. The planform for test case 3. Dashed lines denote the modified 
planform for our calculations. 
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variation of the trail ing-edge height over the x,y-plane, with the leading edge 
of the wing remaining on the x*,y-plane over its entire length. Total twist var- 
iation between the root and the tip is 6 degrees. 

The calculations were made for the upper surface at Mach numbers of 
M =0.5 and 0.99 by using the experimental pressure distributions given in 
[21] as Table III-230 and Table III-28, respectively. The external velocity 
components were obtained by using the procedure described in Appendix A and 
by using the velocity program described in Appendix C. Since the experimental 
pressure distribution does not extend over the entire chord-length, the 
boundary-layer calculations were started at x/c = 0.025 and were continued 

up to (x/c) = 0.975. The transition point was assumed at (x/c) = 0.10. 

6 

The unit Reynolds number in both cases was 4.9 x 10 . The spanwise extent 
of calculation is from y = 0.76 meters to y = 6.6 meters and was again limited 
by the available pressure measurements. For this reason we have extrapolated 
the experimental data to get the data on the root section. For both cases, 
w 0 was negative on the entire wing; as a result calculations were started 
at the tip section and were continued all the way to the root section. Note 
that in terms of the spanwise boundary-layer coordinate z, y = 0.76 meters 
corresponds to z = 0, and y = 6.6 meters corresponds to z - 0.92. 

The computed chordwise and spanwise values of c f and 6* are shown 
for both cases in figures 18 to 23. Note that for = 0.50, the chordwise 
local skin-friction distribution at z = 0.7b is quite different from those 
c^-distributions at z = 0.52 and 0.28. Similarly, the chordwise variation 
of displacement thickness for the same Mach number at each spanwise station 
are quite different from each other. The large variation in both skin- 
friction and displacement thickness distributions in spanwise direction is 
partially due to the experimental pressure distribution selected for this 
sample calculation. That is, the flow was obviously separated on the out- 
board panels of the wing, and it was the separated flow pressure distribu- 
tion which was specified for boundary-layer calculations. For = 0.99, 
the chordwise local skin-friction and the displacement thickness distribu- 
tions for three different spanwise stations are qualitatively similar, as 
one would expect them to be for properly designed wing with attached flow. 


49 




Figure 18. Variation of the local skin-friction in the chordwise direction for M ro = 0.50. (a) z = 0.76. 

The dashed lines denote the laminar flow solutions and the solid lines “denote the turbulent 
flow solutions. 
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Figure 19. (c) z = 0.28 














CONCLUDING REMARKS 


According to the studies presented in this report, the three-dimensional 
boundary layers on swept wings can be computed very effectively by our 
method. The geometry package represents the wing analytically and accounts 
for its geometric parameters. However, the method should be tested further 
and additional capabilities should be added in order to make the method 
more versatile. In particular: 

1. The marching procedure should be further tested. Although the sample 
calculations are very encouraging, additional test cases must be made 
for different pressure distributions on different wings to check the 
numerical calculations for strong negative cross-flow velocity profiles. 

2. The present procedure to calculate the attachment-line flow near the 
wing leading edge should be modified to handle cases where the attach- 
ment streamline deviates considerably from a constant chordwise loca- 
tion along the span. 

3. The boundary-layer method should be coupled to an inviscid code so 
that inviscid flow and viscous flow calculations can be made at the 
same time. 



APPENDIX A 


CALCULATION OF THE EXTERNAL VELOCITY COMPONENTS 
FROM A GIVEN PRESSURE DISTRIBUTION 


Because the experimental pressure distribution is very seldom given with 
sufficient accuracy or in sufficient detail to calculate the inviscid velocity 
components from Euler's equations, we resort to the local sweep theory. The 
local sweep theory is known to give reasonable results when applied to regions 
of high aspect ratio wings that are outside the influence of root and tip 
effects. If the spanwise pressure gradient is weak, as on the midportion of 
a swept wing, the accuracy of the approximate method is almost exact. In 
regions of root and tip influence a correction to the local sweep theory is 
applied in order to account for the effective reduction of the sweep angle. 

Our procedure is explained below. 

Consider the resultant velocity vector u s in the tangent plane at a 
point P on the wing surface (figure Al). The basic assumption for the 
local sweep theory is that projections of both the freestream velocity vector 
u ro and the local velocity vector u s upon the z-axis are equal, or 

Up = u $ sing = u^ sinx (Al ) 

Here x is the local sweep angle obtained from 



From the parallelogram law for the addition of vector components we obtain 

u u. cos g 

u u sine v 


w e _ u s / sing sine — cose cosg 
u u \ sine 

CO oo 

Elimination of g from (A3) and (A4) gives 


(A4) 
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Figure A1 . 



Velocity vector in the tangent plane at a point P 
a is the sweep angle. 

u e _V(u s /uJ 2 - Sin 2 x 
u „ sine 


on the wing. 


(A5) 


w e . u e 

— = sinx — — cose 


The total velocity ratio (u /u^) is calculated from 




P t/ P T 2 (1 + 


(y-1 )/y 


( y - 1)MV2 


(A6) 


(A7 ) 


with Pjj and Pj^ denoting the values of total pressure before and after 
the shock, respectively, and C is the pressure coefficient, C = (P — Pj)/ 
(l/2pu oo ). Equation (A7) is valid for an adiabatic flow through a shock wave, 
but since the total pressure ratio across the shock is seldom known, it is 
set equal to one. This approximation introduces only a few percent error into 
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the velocity calculations because the, total pressure jump across a swept 
shock is small even for freestream Mach numbers approaching unity. The total 
pressure ratio must also remain close to one for the first-order boundary- 
layer theory to be valid in front and behind the shock wave. 


Equations (A5) to (A6) are approximately valid for the root and tip 
regions if the local sweep angle is replaced by an effective sweep angle. 


x e = 


F r,t x r,t 


(A8) 


Here a denotes the local sweep angle, \ . denote the root or tip sweep 

■ 5 ^ 

angle for the given z-coordinate line. F . are the spanwise interpolation 

r jl 

factors for the root and tip, respectively. They are given in graphical form 
in [ TO] as a function of nondimensional spanwise distance in terms of root or 
tip chord. The curve fit for the above interpolation factors are: 



52 

90(y/c r ) + 7 

/ 0.24 

\28(b -y)/c t 



+ 3 


0 <f-< 0. 5 

c r 




0.75 


(A9) 

(A10) 


where b is the semi span, and c r and c^ are the root chord and the tip 
chord, respectively. Since no overlap of (A9) and (A10) is allowed, there is 
a lower limit of wing aspect ratio to which the local sweep theory with cor- 
rection can be applied. For example, for a wing with no taper, the wing aspect 
ratio must be greater than or equal to 2.5. As a consequence of the effective 
reduction of sweep angle, the calculation of velocity components in a nonortho- 
gonal coordinate system must reflect the fact that the angle between the 
coordinate lines is effectively increased. Simple geometrical consideration 
of a sheared wing shows that the effective coordinate angle should be expressed 
as 

sin A 

cos6 e = sTneT cose (A11 > 

To use equations (A5) and (A6) near the wing root or tip A is replaced by 
A e and 0 is replaced by 0 g . 
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APPENDIX B 


CALCULATION OF THE COORDINATE SYSTEM AND ITS GEOMETRIC PARAMETERS 


Coordinate System 

We assume the wing to be defined in the Cartesian coordinate system 
x, y, z. Here x is in the general direction of the airplane longitudinal 
axis, y is in the spanwise direction and z is normal to the xy-plane. 

We shall also assume that the wing definition is given by a number of air- 
foil sections in planes y = const., which involves the specification of 
sets of x is z 1 for constant values of y y The boundary layer is to be 
calculated at a fixed number of nondimensional chordwise locations U/c) 

r 

for a given number of spanwise stations y $ . The coordinate system for 
boundary-layer calculation is defined by the lines (S/c) p = const, and 
y s = const, on the wing surface and the surface normals. Figure Bl shows 
the plan view of a wing with our notation. The choice of the coordinate 
system for boundary- layer calculations is dictated by the fact that aero- 
dynamic data Is usually given in terms of percent of chord and percent of semi 
span location, and by the almost uniform coverage of the wing with the chosen 
coordinate net. 


To find the geometrical properties of the chosen coordinate lines, we 

need relationships between the wing defining system x, y, z and the surface 

imbedded coordinate system (c/c) , y . 

P s 


Considering y = const, planes, the following relationships hold (see 
Figure B2) 


where 


§ = ^ C (x - x £ ) cosa - (z - z % ) Sina] 

r 1 

£ = c [(x - Xj sina + (z - Zj,) cosa] 


(Bl) 

(B2) 


c = -** )2 + (5 t 




(B3) 

(B4) 
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Figure B1 . Plan view of a wing with our notation. 

Here the subscripts i and t refer to the points at the leading edge 
and trailing edges, respectively, and c is the local chord length (maximum 
length line). 

Considering figure B2, it is obvious that ( 5 /c) = const, are accept- 
able coordinate lines as long as we can treat the upper and lower surfaces 
separately. However, if the stagnation point s is located in the lower 
surface, calculations for what we now call the upper surface contain a portion 
of the lower surface and the meaning of coordinate £/c becomes ambiguous. 

To remedy this situation and also to stretch the coordinate near the leading 
edge where flow quantities vary rapidly, we employ the following transforma- 
tion of the independent variable c/ c: 
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Figure B2. Notation for the airfoil section for a given y. 

§ = 1 - cos<$> (B5) 

Here <}> = 0 corresponds to the leading edge, and <j> = tt/ 2 to the trailing 

edge. If the upper surface is being calculated, the points on the lower surface 
correspond to negative values of <f> and vice versa. On the coordinate line 
cj> = const, we use as the independent variable the nondimensional form of y 
given by y = y/b where b is a scaling constant for the convenience of the 
user of the program. If b denotes the length of semi span, y will vary 
from 0 to 1 . The variables y constitute a nonorthogonal coordinate 
system embedded in the wing surface. It should be noted that the boundary- 
layer coordinates in the main text make use of a different notation, namely 
x = <}) and z = y in the present notation. 
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Geometric Parameters 


Once the coordinate system is selected, it is necessary to calculate its 
geometric parameters, such as the metric coefficients, and the geodesic 
curvatures of the coordinate lines. In what follows we designate quantities 
associated with the chordwise variable with subscript one and those associ- 
ated with the spanwise variable with subscript two. 


The metric coefficient along a curve in space is given by 



(B6) 


with P denoting a parameter. Taking P = 4 along the curves y = const., 
we can write (B6) for h-j : 



(B7 ) 


Similarly, taking P = y along the curves <t> = const., we can write 


h 


2 = b 2 + 



(B8) 


The derivatives in (B7) and in (B8), namely (3x/a<|>) , ( 8z/3 4 ) )y » (3x/3y) , 
and (3z/3y)^ are obtained as by-products of spline-fitting the points 
along the chordwise and spanwise directions at the pivotal points. 


The unit tangent vector ? along a curve is given by 


The unit tangent vector 


t - dr _ dr 1 1 dr 

1 = ds - W Xds/dPT = FdP 

tl along the curve y = const. 



is 


(B9) 


(BIO) 


where i, j, k are unit vectors in the coordinate directions x, y, 5, 
respectively. 
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The unit tangent vector % along the curve <p = const, is 


t = J— 
z 2 Ti^ 


.(»), ’ * « * («)/ 
The angle between the coordinate lines is then 


(BID 


cose = ti • tg = 


(si#) * (#)(#: 

The curvature of a curve in space is given by 


/(h-|h 2 ) 


t, _ d? _ dt 1 1 d? 

K is W (ds/dP) h dP 


(B12) 


(B13) 


The geodesic or tangential curvature Kg of a curve on the surface can be 
obtained from 


K g = (? * n) • K 

Here n is the vector normal to the surface which by definition is 

sinen = t-| x t 2 


or 


-> = 1 
n h-j h 2 s' i n e 


-b /lih - 

/ 3 X 

3 z 

dx 

liV- + Ja 

a k l 

_ b UJ 

\» 4 > 

ay 

ay 

3 \3 

*rj 


With the use of (B9), the expression ( B 1 3 ) can be written as 


(B14) 


(B15) 


(B16) 


t 1 d 2 r 1 dr dh _ 1 / d 2 x . d 2 y . 

In _ ~~yy — ^ tTFr tft - — ?r I — nr l + — n .1 


, + d 2 z ; 

7^ h 3 dP ^'7U 2 ' " dP 2 J dP 2 / 


1_ /dx i . dy , dz \( d 2 x dx d^£ d£ . d 2 z dz' 

hA* 7 + + 3F^dP + ^# + ^dF, 


The geodesic curvature Kg^ for curve y = const, is 


Kg-j = -(t-j x n) • t-j 


(B17) 


(B18) 


The minus sign on the right-hand side of ( BIB) is introduced to obtain 
Kg i = ~(l/h-| h 2 ) (ah-|/ 3 y) in the case of orthogonal coordinate system. With 
<j> as the parameter, the expression for K-| is 
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K i M? 


[_3<f> 


2 - 2 - 

4- 1 + 14 k 


3 (j) 


1 

* 


3X ^ ^ 3 Z j I o a o a a a o^. 


2- - 2- - N 
3 X 3x . 3 Z 3Z 


3 tf> 


3 c j > 


8*- 8 * ^7 


(B19) 


Introducing the expressions (BIO), (B16) and (B19) into (B18) gives, after 
simplifications 


K, 
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-1 


( - _ _ _\ / 2- - 2- 

3X 3Z 3x 3 Z \ / 3 X 3Z 3 Z 3X 


3* a+ 2 3 4> 


h|h 2 sine v* 3y 3y a V\3« 

The geodesic curvature for a curve $ = const, is given by 

Kg 2 = ^2 * ^ * ^2 

With y as the parameter, the expression for IC, is 
K = - 1 - 

2 h 2 
ho 


r 2 - ?- 

l 


r ?- - 

2- -1 

i + k 

_3y sy . 

3 X . . . . , 3Z . 

3y 1 + bj + W k 

| 3 X 3X 3 Z 3Z 

" h T 

h 2 

Lay 2 

3y 2 3y _ 


(B20) 


(B21) 


(B22) 


The expression for the geodesic curvature Kg 2 is obtained by substitution 
of (BIT), (B16) and (B22) into (B21): 
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h-|h 2 sine 


3x 3 z 
3tj> 3y 


3x 3 Z 
3y 3<j> 



3 2 Z 3x\ , ,2 / 3 2 X 3x , 3 2 z 3 z 

b W*' 77^ 


(B23) 

As discussed in the beginning of this section, the program reads in n 
pairs (n can be variable along span) of x^ , z^ values at each spanwise 
station y*. Then the use of equations (Bl) and (B5) gives the corresponding 

J 

cj>.. Since the desired chordwise outputs are read in as (c/c) , we again 

1 ^ r ^ 

use equation (B5) to obtain $ . Next, the tables of x.. and z^ vs <j> - 
at each y- are fitted with cubic splines to obtain interpolated values of 

J 

x , z , and 3z /3cj> corresponding to <|> . Here the cubic splines are not 

r r r r 

used to obtain the second derivatives because of inherent inaccuracies. 

Instead, we use a Fourier fit with sigma-smoothing. The end result is 
— 2-2 

smoothed 3Zp/3<|> values and 3 Zp/3<j> . Because of relationships (Bl ) 
and (B5), the derivatives of x can be calculated as follows: 
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tana 


3z rt c sin<f> 

— B. + IP 

9<j> COSa 



tana 


2 - 

3 Z p C COScjjp 

7J + COSa 

d < P 


(B24) 


(B25) 


To obtain the derivatives with respect to y at desired output locations 
y s , we spline-fit the Xp and Zp values at constant cf>p vs y^.. The 
first derivatives are again not used directly, instead they are put through 
the Fourier smoothing procedure and simultaneously interpolated for values 
at y s together with the second derivatives. To obtain the x p and Zp 
values and their derivatives with respect to 4> at y s cubic spanwise 
interpolation is used. 

After the derivatives are determined the quantities h-j , hg, <j>, Kg^ 
and Kg 2 are calculated from the relationships given previously. To cal- 
culate the remaining quantities K^ and requires extra care at the 

leading edge because the terms K-j and 1/h-j 9 e/3 (f are large and of 
opposite sign. In effect, the sum of these two terms is similar to a third- 
order derivative. The best results were obtained by sigma-smoothing the 
calculated K-j and values and applying the same technique to the cal- 
culation of e-derivatives. It is also helpful to extend the wing defini- 
tion input data a long way around the leading edge on the opposite surface 
if boundary-layer calculations involve the leading edge. 
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